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ALGEBRAIC SYSTEMS BIOLOGY: 

A CASE STUDY FOR THE WNT PATHWAY 

ELIZABETH GROSS, HEATHER A. HARRINGTON, ZVI ROSEN, AND BERND STURMFELS 


Abstract. Steady state analysis of dynamical systems for biological networks give rise to 
algebraic varieties in high-dimensional spaces whose study is of interest in their own right. 
We demonstrate this for the shuttle model of the Wnt signaling pathway. Here the variety 
is described by a polynomial system in 19 unknowns and 36 parameters. Current methods 
from computational algebraic geometry and combinatorics are applied to analyze this model. 


1. Introduction 


The theory of biochemical reaction networks is fundamental for systems biology 13, 27 


It is based on a wide range of mathematical fields, including dynamical systems, numerical 
analysis, optimization, combinatorics, probability, and, last but not least, algebraic geometry. 
There are numerous articles that use algebraic geometry in the study of biochemical reaction 
networks, especially those arising from mass action kinetics. A tiny selection is @01 12,22,25 


We here perform a detailed analysis of one specific system, namely the shuttle model for the 


Wnt signaling pathway, introduced recently by MacLean, Rosen, Byrne, and Harrington 17 


Our aim is twofold: to demonstrate how biology can lead to interesting questions in algebraic 
geometry and to apply state-of-the-art techniques from computational algebra to biology. 


The dynamical system we study consists of the following 19 ordinary differential equations. 
Their derivation and the relevant background from biology will be presented in Section [2} 


X\ 

x-2 

X'3 

X 4 

x 5 

x 6 

x 7 

fll "^8 ^T6 

^ ' Xg = ~X 17 

Xio 

Xu 

Xl2 = ~Xl3 
X U 
Xl5 
Xl8 
Xl9 


~k iTi + k 2 X 2 

k\X\ — (k 2 + k 2 ^)x 2 + k 27 X3 — k3X 2 x A + (k 4 + k^)xi A 
k 2 §X 2 — k 27 X3 — k\4X3X5 + (&15 + ^ 16)^15 
-k 3 x 2 X4 - k 9 x 4 x w + k A x 14 + k 8 x 16 + (k 10 + k n )x 18 
-k 2S x 5 + k 29 x 7 - k 6 x 5 x 8 + k 5 x u + k 7 x 16 
-k lA X3XQ - k 20 XQXu + fci 5 xi5 + ki 9 x l7 + (k 2 1 + k 22 )x 19 
k 2S x 5 - k 29 x 7 - k 17 x 7 x 9 + k 16 x l5 + k 18 x 17 
-k e x 5 x 8 + (hr + k 8 )x w 
— k\ 7 x 7 x 9 + (kis + ki 9 )xi 7 

k \2 — (kl 3 + ^ 30)^10 — kgX 4 Xio + k 8 \X 11 + kioXis 

— k 2 3Xn + ^ 30^10 — ^ 3 lUl — k 29 XQXu — k 2 4;Xi\Xi 2 + k 2 §X\3 + ^ 21^19 

k 2 4X\\X\ 2 + k 2 $X\3 
k 3 x 2 x 4 - (k A + k 5 )x 14 
ki4X 3 x 6 - (k 15 + k w )x 15 
k 9 x 4 x 10 - (ho + kn)xis 
k2oX 6 Xn - (k 2 1 + k 22 )x 19 
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The quantity x\ is a differentiable function of an unknown t, representing time, and Xi(t) is 
the derivative of that function. This dynamical system has five linear conservation laws: 


( 2 ) 


0 = (xi + X 2 + x 3 + X 14 + X15) - Cl 

0 = (aq + x§ + Xq + Xj + £14 + £15 + Xi6 + X17 + + aqg) — C2 

0 = (x 8 + x 16 ) - c 3 

0 = (xg + X17) - C 4 

0 = {x 12 + ^13) — c 5 


The 31 quantities ki are the rate constants of the chemical reactions, and the five c* are the 
conserved quantities. Both of these are regarded as parameters, so we have 36 parameters 
in total. Our object of interest is the steady state variety , which is the common zero set of 
the right hand sides of ([!]) and (|2|. This variety lives in A' 19 , where K is an algebraically 
closed field that contains the rational numbers Q as well as the 36 parameters ki and c % . If 
these parameters are fixed to be particular real numbers then we can take K = C, the field 
of complex numbers. If it is preferable to regard k = (Aq,..., fc 31 ) and c = (ci,..., C 5 ) as 
vectors of unknowns, then K = Q(k, c) is the algebraic closure of the rational function field. 
In this latter setting, when all parameters are generic, we shall derive the following result: 

Theorem 1.1. The polynomials in |?])-|^|) have 9 distinct zeros in K 19 when K = Q(k, c). 


By analyzing the steady state variety, we can better understand the model, which is non¬ 
linear, and thus the biological system. The aim is to predict the system’s behavior, offer 
biological insight, and determine what data are required to verify or reject the model. Here 
is a list of questions one might ask about our model from the perspective of systems biology. 


Biological Problems. These are labeled according to the section that will address them. 

4. For what real positive rate parameters and conserved quantities does the system exhibit 
multistationarity? This question is commonly asked when using a dynamical system 
for modeling a real-world phenomenon. When modeling a process that experimentally 
appears to have more than one stable equilibrium, multistationary models are preferred. 

5. Suppose we can measure only a subset of the species concentrations. Which subsets can 
lead to model rejection? If all species are measurable at steady state, then we can substi¬ 
tute data into the system ([Tj) , and check that all expressions iq are close to zero. If only 
some Xi are known, we still want to be able to evaluate models with the available data. 

6 . Give a complete description of the stoichiometric compatibility classes for the chemical 
reaction network. A stoichiometric compatibility class is the set of all points accessible 
from a given state via the reactions in the system. This question relates more closely to 
the dynamics of the system, but also has ramifications for the set of all steady states. 

7. What information does species concentration data give us for parameter estimation? 
In particular, are the parameters identifiable? Identihability means that having many 
measurements of the concentrations x can determine the reaction rate constants k. If not 
identifiable, we will explore algebraic constraints imposed by the species concentration 
data. This question is relevant for complete and partial steady-state data (usually noisy). 


These questions are open challenges for medium to large models in systems biology and 
medicine 13,27. The book chapter 16 illustrates standard mathematical and statistical 
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methods for addressing these questions, with Wnt signaling as a case study. Here, we examine 
these questions from the perspective of algebraic geometry. The aim is to provide insight into 
global behavior by applying tools from nonlinear algebra to synthetic and systems biology. 
Below are the algebraic problems underlying the four biological problems listed above. 


Algebraic Problems. 

4. Describe the set of points (k, c) 6 x U>o such that the polynomials 0-0 have two 
or more positive zeros x e M^g. When is there only one? Identify the discriminant. 

5. Which projections of the variety defined by (Jl| into coordinate subspaces of K 19 are 
surjective? Equivalently, describe the algebraic matroid on the ground set (xi,... , X 19 }. 

6 . The conservation relations (J2]) specify a linear map x '■ ®1 19 —> ®1 5 , x4c. Describe all the 
convex polyhedra X -1 (c) D M > 9 0 where c runs over the points in the open orthant ®>o- 

7. a. Complete data: Describe the matroid on the ground set {k\ } k- 2 ,..., k 3 i} that is defined 

by the linear forms on the right hand sides of (jl|, for fixed steady-state concentrations. 

b. Partial steady-state data without noise: Repeat the analysis after eliminating some of 
the x-coordinates. 


c. Partial steady-state data with noise: For the remaining x-coordinates, suppose that we 
have data which are approximately on the projected steady state variety. Determine a 
parameter vector (k, c) that best fits the data. 


In this paper we shall address these questions, and several related ones, after explaining the 
various ingredients. A particular focus is the exchange between the algebraic formulation 
and its biological counterpart. Our presentation is organized as follows. 

In Section [2] we review the basics on the Wnt signaling pathway, we recall the shuttle model 
of MacLean et al. [17], and we derive the dynamical system ((TJ) (|2]) . In Section[3]we establish 
Theorem 1.1 and we examine the set of all steady states. This is here regarded as a complex 
algebraic variety in an affine space of dimension 55 = 19 + 31 + 5 with coordinates (x, k, c). 

In Sections [4] [5j [6j and [ 7 ] we address the four problems stated above. The numbers of 
the problems refer to the respective sections. Each section starts out with an explanation 
of how the biological problem and the algebraic problem are related. The rationale behind 
Section [4] is likely to be familiar to most of our readers, given that multistationarity has been 
discussed widely in the literature; see e.g. (4, 221. On the other hand, in Section [5] we employ 
the language of matroid theory. This may be unfamiliar to many readers, especially when 
it comes to the algebraic matroid associated with an irreducible algebraic variety. Section [6] 
characterizes the polyhedral geometry encoded in the conservation relations (|2]). This is a 
case study in the spirit of |25, Figure 1]. Section [ 7 ] addresses the problems of parameter 
identihability and parameter estimation. Finally, in Section [8] we return to the biology, and 
we discuss what our findings might imply for the study of Wnt signaling and other systems. 
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2 . From Biology to Algebra 


Cellular decisions such as cell division, specialization and cell death are governed by a rich 
repertoire of complex signals that are produced by other cells and/or stimuli. In order for 
a cell to come to an appropriate decision, it must sense its external environment, communi¬ 
cate this information to the nucleus, and respond by regulating genes and producing relevant 
proteins. Signaling molecules called ligands, external to the cell, can bind to proteins called 
receptors, initializing the propagation of information within the cell by molecular interactions 
and modifications (e.g. phosphorylation). This signal may be relayed from the cytoplasm 
into the nucleus via molecules and the cell responds by activation or deactivation of gene(s) 
that control, for example, cell fate. The complex interplay of molecules involved in this in¬ 
formation transmission is called a signaling transduction pathway. Although many signaling 
pathways have been defined biochemically, much is still not understood about them or how a 
signal results in a particular cellular response. Mathematical models constructed at different 
scales of molecular complexity may help unravel the central mechanisms that govern cellular 
decisions, and their analysis may inform and guide testable hypotheses and therapies. 

In this paper, we focus on the canonical Wnt signaling pathway, which is involved in cel¬ 
lular processes, both during development and in adult tissues. This includes stem cells. 
Dysfunction of this pathway has been linked to neurodegenerative diseases and cancer. Con¬ 
sequently, Wnt signaling has been widely studied in various organisms, including amphibians 
and mammals. Researchers are interested in how the extracellular ligand Wnt affects the 
protein /3-catenin, which plays a pivotal role in turning genes on and off in the nucleus. 

The molecular interactions within the Wnt signaling pathway are not yet fully understood. 
This has led to the development and analysis of many mathematical models. The Wnt 
shuttle model jlT| includes an abstraction of the signal transduction pathway (via activa¬ 
tion/inactivation of molecules) described above. The model also takes into account molecules 
that exist, interact and move between different compartments in the cell (e.g., cytoplasm 
and nucleus). Biologists understand the Wnt system as either Wnt off or Wnt on. How¬ 
ever, such a scenario is rarely binary (i.e., different concentration levels of Wnt may exist) 
and inherently depends on spatial movement of molecules. The Wnt shuttle model includes 
complex interactions with nonlinearities arising in the equations. In particular, it includes 
both the Wnt off and Wnt on scenarios, by adjusting initial conditions or parameter values. 
The biology needed to understand the model can be described as follows. See also Table [TJ 

Wnt off: When cells do not sense the extracellular ligand Wnt, /3-catenin is degraded (bro¬ 
ken down). The degradation of /3-catenin is partially dependent on a group of molecules 
(Axin, APC and GSK-3) that form the destruction complex. Crucially, the break down of 
/3-catenin occurs when the destruction complex is in an active state; modification to the 
destruction complex by proteins, called phosphatases, changes it from inactive to active. 
Additionally, /3-catenin can degrade independent of the destruction complex. Synthesis of 
/3-catenin occurs at a constant rate. 

Wnt on: When receptors on the surface of a cell bind to Wnt, the Wnt signaling transduc¬ 
tion pathway is initiated. This enables /3-catenin to move into the nucleus where it binds 
with transcription factors that regulate genes. This signal propagation is mediated by the 
following molecular interactions. After Wnt stimulus, the protein Dishevelled is activated 
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near the membrane. This in turn inactivates the destruction complex, thereby preventing 
the destruction of /3-catenin, allowing it to accumulate in the cytoplasm through natural 
synthesis. Throughout the molecular interactions in the signaling pathway, intermediate 
complexes can form (e.g., /3-catenin bound with Dishevelled). 


Space: The location of molecules plays a pivotal role: /3-catenin moves between the cytoplasm 
and the nucleus (to reach target genes and regulate them). Dishevelled and molecules that 
form the destruction complex shuttle between the nucleus and the cytoplasm. However, it 
is assumed that only the inactive destruction complex can shuttle (since in the cytoplasm it 
would be bound to /3-catenin). Phosphatases exist in both the nucleus and the cytoplasm 
but the movement across compartments is not included in the model. Symmetry of reactions 
is assumed if the species exist in both compartments. Intermediate complexes are assumed 
to be short-lived, or not large enough for movement across compartments. 


The Wnt shuttle model of 17 has 19 species whose interactions can be framed as biochemical 
reactions. These species correspond to variables x\,...,xig in our dynamical system ( |I|). 
Namely, x % represents the concentration of the species that is listed in the it\i row in Table]!] 


Variable 

Species 

Symbol 


Dishevelled 

D 

X\ 

Dishevelled in cytoplasm (inactive) 

Di 

X 2 

Dishevelled in cytoplasm (active) 

D a 

X 3 

Dishevelled in nucleus (active) 

D a n 


Destruction complex (APC/Axin/GSK3/3) 

Y 

X 4 

Destruction complex in cytoplasm (active) 

Y a 

x 5 

Destruction complex in cytoplasm (inactive) 

Yi 

x 6 

Destruction complex in nucleus (active) 

Y an 

X7 

Destruction complex in nucleus (inactive) 

Y 

1 m 


Phosphatase 

p 

x 8 

Phosphatase in cytoplasm 

p 

Xg 

Phosphatase in nucleus 

Pn 


(3— catenin 

X 

XlQ 

/3-catenin in cytoplasm 

X 

Xu 

/3-catenin in nucleus 

X n 


Transcription Factor 

T 

Xl 2 

TCF (gene transcription in nucleus) 

T 


Intermediate complex 

C 

Xl3 

Transcription complex, /3-catenin: TCF in nucleus 

C X T 

Xu 

Intermediate complex, /3-catenin: dishevelled in cytoplasm 

Cyd 

Xl5 

Intermediate complex, destruction complex: dishevelled in nucleus 

CyDu 

X W 

Intermediate complex, destruction complex: phosphatase in cytoplasm 

Cyp 

xn 

Intermediate complex, destruction complex: phosphatase in nucleus 

CyPu 

XlS 

Intermediate complex, /3-catenin: destruction complex in cytoplasm 

C x y 

Xig 

Intermediate complex, /3-catenin: destruction complex in nucleus 

DxYn 


Table 1: The 19 species in the Wnt shuttle model. 
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The second column in Table [T| indicates the biological meaning of the 19 species. The symbols 
in the last column are those used in the presentation of the Wnt shuttle model in [l7 |. 

The 19 species in the model interact according to the 31 reactions given in Table [2} Each 
reaction comes with a rate constant fcj. These are the coordinates of our parameter vector k. 


Reaction 

Explanation 

ki 

Xi x 2 

ki 

(In)activation of dishevelled, depends on Wnt 

k 3 , k 5 

X 2 + X 4 x 14 - > X 2 + X 5 

Destruction complex active —¥ inactive 

fcg 

X5 + X 8 <=> Xi 6 ->■ X 4 + X 8 

k 7 

Destruction complex inactive —>■ active 

kf) s fcl1 , rt 

x 4 + Xio Xis - > x 4 T v 

kio 

Destruction complex-dependent /Tcatcnin degradation 

rt fcl2 

0 - ¥ X\q 

/3-catenin production 

k 13 ft 

Xio -^0 

Destruction complex-independent /3-catenin degradation 

ki4 kie 

X 3 + Xq X 15 -!• X 3 + X 7 

k 15 

Destruction complex active —>■ inactive (nucleus) 

k i7 k 19 

X 7 + Xg X 47 - ¥ Xq + Xg 

kis 

Destruction complex inactive —¥ active (nucleus) 

fc 2° k 22 , ft 

Xq + Xu 5==^ XiQ - ¥ X 6 + (/) 

k 2 \ 

Destruction complex-dependent /3-catenin degradation (nucleus) 

k 23 ft 

Xu -> 0 

Destruction complex-independent /3-catenin degradation (nucleus) 

k 2 4 

Xu + X 12 1 =^ x 13 

k 2 5 

/3-catenin binding to TCF (nucleus) 

k 2 Q 

X 2 l=tx 3 

k 27 

Shuttling of active dishevelled 

k 2 g 

x 5 X 7 

k 29 

Shuttling of inactive-form destruction complex 

k30 

Xio Xu 

ksi 

Shuttling of /3-catenin 

Table 2: 

The 31 reactions in the Wnt shuttle model. 


The 31 reactions in Table [2] translate into a dynamical system x = T(x: k). Here T is a 
vector-valued function of the vectors of species concentrations x and rate constants k. The 
choice of T is up to the modeler. In this paper, we assume that T represents the law of 
mass action 13, §2.1.1], This is precisely what is used in for the Wnt shuttle model. 
The resulting dynamical system is ([!]). We refer to 4,7,12,22,25 and their many references 


for mass action kinetics and its variants. In summary, Table |2 translates into the dynamical 
system Q under the law of mass action. The five relations in ([2]) constitute a basis for the 
linear space of conservation relations of the model in Table [2] assuming mass action kinetics. 
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We refer to xi,...,x± 9 as the species concentrations , ki ,..., fc 3 i as the rate parameters , 
and Ci,... ,c 5 as the conserved quantities. We write x, k and c for the vectors with these 
coordinates. As is customary in algebraic geometry, we take the coordinates in the complex 
numbers C, or possibly in some other algebraically closed field K containing the rationals Q. 

Our aim is to understand the relationships between x, k and c in the Wnt shuttle model. To 
this end, we introduce the steady state variety S C C 55 . This is the set of all points (x, k, c) 
that satisfy the equations — ... — ± 19 = 0 in 0 along with the five conservation laws 
in (|2]). We write our ambient affine space as C 55 = C^ 9 x C^ 1 x Cjh This emphasizes the 
distinction between the species concentrations, rate parameters, and conserved quantities. 

3. Ideals, Varieties, and Nine Points 

We write / for the ideal in the polynomial ring Q[x, k] = Q[aq,... x\ 9 , ki,... k 8 i\ that is 
generated by the 19 polynomials ±j on the right hand side of ([Tj). Five of these generators 
are redundant. Indeed, the conservation relations (|2]) give the following identities modulo /: 

±1 + ±2 + ±3 + ±14 + ±15 = ±8 + ±16 = ±9 + ±17 = ±12 + ±13 = 

±4 + ±5 + ±6 + ±7 + ±14 + ±15 + ±16 + ±17 + ±18 + ±19 = 0. 

For instance, the polynomials ± 13 , ± 15 , ± 16 , ±17 and ±19 are redundant because they can be 
expressed as negated sums of other generators of I. Hence / is generated by 14 polynomials. 
The variety V(I) lives in the 50-dimensional affine space C^ 9 x C^ 1 , and it is isomorphic to 
the steady state variety S C C 55 . A direct computation using the computer algebra package 
Macaulay2 |TT] shows that V(I) has dimension 36. Hence the affine ideal / is a complete 
intersection in Q[x, k]. Furthermore, using Macaulay2 we can verify the following lemma. 

Lemma 3.1. The ideal I admits the non-trivial decomposition I = J m DJ e; where I e = I : (aq) 
and I m — I + (aq), both of these components have codimension 14, and I e is a prime ideal. 

The ideal I m is called the main component, while I e is called the extinction component , since 
it reflects those steady states where a number of the reactants “run out.” Both of these ideals 
live in Q[x, k], and we now present explicit generators. The extinction component equals 

Ie — { x \,X' 2 , £ 3 , £ 5 , X 7 , X 14 , X 15 , X\q, X\ 7 , k 89 X\Q — ( k 23 + & 3 l)xn — k 22 X\ 9 , 

^13^10 + k 28 Xn + kuXis + k 22 X\ 9 — &12, k 2 4XuX\ 2 — ^ 25 X 43 , 

k 2 ox e xii - (k 2 1 + k 22 )x 19 , k 9 x A x w - (Aq 0 + Aqi)aq 8 ). 

The ideal I e is found to be prime in Q[x, k]. The main component equals 

Im = (kl6Xi5 — kigXn, k^Xu — k 8 X 16, ^30^10 — (&23 + ^3l)^ll — k 22 X\ 9 , 

ki^Xio + k 28 Xn + kiiX\ 8 + k 22 xi 9 — k 72 , k 28 x§ — k 29 x 7 , k 26 x 2 — k 27 x 8 , 
k\X\ - k 2 x 2 , k 2A XuX l2 - k 25 x 18 , k 20 x 6 x u - (k 21 + k 22 )x 19 , 
kgX 4 X 10 - (Alio + kn)x 18 , k 17 x 7 x 9 - (k 18 + k 19 )x 17 , k 6 x 5 x 8 - (k 7 + k 8 )x i 6 , 
kux 8 XQ — ki 5 xia — k\ 9 Xi 7l k 8 x 2 X 4 — Aqaq 4 — k 8 x 7 §, 

{k4k 9 k 8 ki4ki§ki 8 k 29 k 29 + kz,kQk 8 k\4ki§k\ 8 k 2 Qk 29 -\- 
k4k 9 k 8 ki4ki§ki 9 k 29 k 29 + k 7) k§k 8 k\4k\§k\ 9 k 2 sk 29 )k\x§x 8 
— (k 8 k^k 7 ki5ki 7 ki 9 k 27 k 28 + k 8 k5k 8 ki5ki 7 ki 9 k 27 k 28 
-\-k 8 k§k 7 kiQki 7 ki 9 k 27 k 28 + k 8 k^k 8 k\Qki 7 k\ 9 k 27 k 28 ) k\X4X 9 ). 
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This ideal is not prime in Q[x, k]. For instance, the variable k\ is a zerodivisor modulo I m , 
as seen from the last generator. Removing the factor k\ from the last generator yields the 
quotient ideal I m : (ki). However, even that ideal still has several associated primes. All of 
these prime ideals, except for one, contain some of the rate constants fcj. 

That special component is characterized in the following proposition. Given any ideal J C 
Q[x, k], we write J = Q(k)[x]J for its extension to the polynomial ring Q(k)[x] in the 
unknowns xi,, Xig over the field of rational functions in the parameters k\,..., kg\. 

Proposition 3.2. The ideal J m = I m D Q[x, k] is prime. Its irreducible variety V ( J m ) C C 50 
has dimension 36; it is the unique component ofV(I m ) that maps dominantly onto C^ 1 . 


Proof. The ideal I m has the same generators as I m but now regarded as polynomials in x with 
coefficients in Q(k). Symbolic computation in the ring Q(k)[x] reveals that I m is a prime 
ideal. This implies that .J m is a prime ideal in Q[x, k], and hence V(J m ) is irreducible. The 
dimension statement follows from the result of Lemma |3.1| that I m is a complete intersection. 
This ensures that V (. I m ) has no lower-dimensional components, by Krull’s Principal Ideal 
Theorem. Finally, V(J m ) maps dominantly onto C^ 1 because J m fl Q[k] = {0}. □ 

Corollary 3.3. The ideal I is radical, and it is the intersection of two primes in Q(k)[x]; 

(3) I = Ie n T m . 


Proof. This follows directly from Proposition 3.2 and the primality of I e in Lemma 3.1 □ 


The decomposition has the following geometric interpretation. We now work over the field 
K = Q(k). All rate constants are taken to be generic. Then V(I) is the 5-dimensional 
variety of all steady states in A' 19 . This variety is the union of two irreducible components, 


V(I) = V(I e ) U V(I m ), 


where each component is 5-dimensional. The first component lies inside the 10-dimensional 
coordinate subspace V(x\, x 2 , x 3 , x 5 , x?, xi 4 , x i5 , x^q, Xn). Hence it is disjoint from the hy¬ 
perplane defined by the first conservation relation x\ + X2 + X3 + Xu + x\$ = C\ . In other 
words, V(I e ) is mapped into a coordinate hyperplane under the map % : K 19 — y A" 5 ,x H- c. 


On the other hand, the second component V (J m ) maps dominantly onto K 5 under y. The¬ 


orem 


( 4 ) 


1.1 states that the generic fiber of this map consists of 9 reduced points. Equivalently, 

x _1 ( c ) n v (7) = x ~\c)nv{Z) 


is a set of nine points in A' 19 . We are now prepared to argue that this is indeed the case. 


Computational Proof of Theorem 0 We consider the ideal of the variety Q in the polyno¬ 
mial ring Q(k, c)[x]. This polynomial ring has 19 variables, and all 36 parameters are now 
scalars in the coefficient field. This ideal is generated by the right hand sides of ([!]) and ([ 2 ]). 
Performing a Grobner basis computation in this polynomial ring verifies that our ideal is 
zero-dimensional and has length 9. Hence Q is a reduced affine scheme of length 9 in A' 19 . 

Fast numerical verification of this result is obtained by replacing the coordinates of k and c 
with generic (random rational) values. In Macaulay2 one finds, with probability 1, that the 






resulting ideals in Q[x] are radical of length 9. We also verified this result via numerical 
algebraic geometry , using the two software packages Bert ini [l] and PHCpack 1261. □ 


4. Multistationarity and its Discriminant 


This section centers around Question 4 from the Introduction: For what real positive rate 
parameters and conserved quantities does the system exhibit multistationarity? This is com¬ 
monly asked about biochemical reaction networks and about dynamical systems in general. 

Mathematically, this is a problem of real algebraic geometry. Writing S for the steady state 
variety in C 55 , we are interested in the fibers of the map 7r k-c : S n M> 0 -> K> o k x M> 0c . 
According to Theorem 0 the general fiber consists of 9 complex points x G C* 9 , when the 
map 7r k ,c is taken over C. But here we take it over the reals M or over the positive reals M >0 . 

In our application to biology, we only care about concentration vectors x whose coordinates 
are real and positive. Thus we wish to stratify k x ®>0,c according to the cardinality of 

(5) 7r k *(k, c) = { (x, k 7 , c 7 ) G S fl : k 7 = k and c 7 = cj. 

This stratification comes from a decomposition of the 36-dimensional orthant k x ®>0,c 

into connected open semialgebraic subsets. The walls in this decomposition are given by the 
discriminant A, a giant polynomial in the 36 unknowns (k, c) that is to be defined later. 

We begin with the following result on what is possible with regard to real positive solutions. 

Theorem 4.1. Consider the polynomial system in 0-0 where all parameters ki and Cj 
are positive real numbers. The set 0 of positive real solutions can have 1,2, or 3 elements. 


Proof. For random choices of (k, c) = (fci,..., fc 31 , Cj,. . ■, c 5 ) in the orthant M> 6 0 , our polyno¬ 
mial system has 9 complex solutions, by Theorem 1.1 For the following two special choices 
of the 36 parameter values, all 9 solutions are real. First, take (k, c) to be the vector 


(1.7182818, 53.2659, 3.4134082,0.61409879,0.61409879, 3.4134082,0.98168436, 0.98168436, 
92.331732, 0.86466471, 79.9512906, 97.932525,1, 3.2654672,0.61699064, 0.61699064, 
37.913879,0.86466471, 0.86466471, 4.7267833, 0.17182818, 0.68292191,1,0.55950727, 
1.0117639,1.7182818,1.7182818, 0.99326205, 0.99326205, 5.9744464,1, 4.9951026, 
16.4733784,1.6006340000000001,1.2089126, 2.7756596399999998). 


The resulting system has three positive solutions x G Mi, 9 0 . Next, let (k 7 ,c 7 ) be the vector 

(0.948166, 7.45086, 5.72974,3.96947, 7.21145, 7.8761,1.87614, 8.11372, 6.21862, 5.24801, 
3.10707,1.08146, 5.22133, 5.84158, .911392,4.28788,4.81201, 9.67849,1.34452, 7.38597, 
6.64451, 7.10229, 8.57942, 5.79076,6.33244,1.53916,1.39658, 0.81673,5.8434, 3.86223, 
7.22696, 1.45438,3.36482,6.06453,4.82045, 3.6014). 

Here, one solution to our system is positive. By connecting the two parameter points above 
with a general curve in M> 6 0 , and by examining in-between points (k 77 , c"), we can construct a 
system with two positive solutions. All computations were carried out using Bert ini [l]. □ 

Remark 4.2. At present, we do not know whether the number of real positive solutions can 
be larger than three. We suspect that this is impossible, but we currently cannot prove it. 

9 




The difficulty lies in the fact that the stratification of M{ 6 0 is extremely complicated. In 
computer algebra, the derivation of such stratifications is known as the problem of real root 
classification. For a sample of recent studies in this direction see |3,6,231. Real root classi¬ 
fication is challenging even when the number of parameters is 3 or 4; clearly, 36 parameters 
is out of the question. The stratification of K* by behavior of ([5j has way too many cells. 


While symbolic techniques for real root classification are infeasible for our system, we can use 
numerical algebraic geometry |9| to gain insight into the stratification of R{ 6 0 . Coefficient- 
parameter homotopies [0j can solve the steady state polynomial system for multiple 

choices of (k, c) quickly. For our computations we use Bertini.m2. This is the Bertini 
interface for Macaulay2, as described in |2|. Each system has 19 equations in 19 unknowns 
and, for random (k, c), each system has 9 complex solutions. Such a system can be solved 
in less than one second using the bertiniParameterHoraotopy function from Bertini.m2. 


Below we describe the following experiment. We sample 10, 000 parameter vectors (k, c) 
from two different probability distributions on R> 6 0 . In each case we report the observed 
frequencies for the number of real solutions and number of positive solutions. We then follow 
these experiments with a specialized sampling scheme for testing numerical robustness. 

Uniform sampling scheme: Here we choose (k, c) uniformly from the cube (0.0,100.0) 36 . 
Sampling 10,000 parameter vectors from this scheme and solving the steady state system for 
each of these parameter vectors in Bertini, we obtained 9, 992 solutions sets that contained 
9 complex points. Solution sets with less than 9 points occur when some paths in the 
coefficient-parameter homotopy fail. We call solution sets with 9 solutions good. 

Integer sampling scheme: Here we select (k, c) uniformly from {1,2,3} 36 . Sampling 10,000 
parameter vectors according to this scheme and solving the corresponding steady state system 
returned 9, 963 good solution sets. Below is a table that records how many of the good 
solution sets had 9, 7, 5, 3 real solutions; all solution sets had 1 positive real solution. 


jf of real solutions 

9 

7 

5 

3 

Freq. for Uniform Sampling 

5,760 

3,675 

544 

13 

Freq. for Integer Sampling 

2,138 

5,181 

2,522 

122 

Table 3: Frequencies for 

die sampling schemes. 


These computations indicate that for most parameter vectors in (0,100) 36 we will see only 
one positive solution to the steady state system. But while the set of parameter vectors that 
result in multiple steady states is not very large, we can give evidence that multistationarity 
is preserved under small perturbations. This is our next point. 


Testing Robustness: Let (k*,c*) be the first point in the proof of Theorem 4.1 For each 


index i G {1,..., 19} we choose yi uniformly from (—0.03 • k *, 0.03 • k*) then set /q = k* Tyt- 
We ran the same process for the c,. Sampling 10,000 parameter vectors this way and solving 
the corresponding steady state systems returned 10, 000 good solution sets, as follows: 


jf of real solutions 

Freq. 


jf of pos. solutions 

Freq. 

9 

9,879 


3 

9,879 

7 

121 


1 

121 


Table 4: Frequencies for testing robustness scheme. 
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In the remainder of this section, we properly define the discriminant A that separates the 
various strata in M> 6 0 . Let A int denote the Zariski closure in C^ 1 x of all parameter vectors 
(k, c) for which (JT|) ([2]) does not have 9 isolated complex solutions and there are no solutions 
with xy = 0 for some i. It can be shown that A int is a hypersurface that is defined over Q, so 
it is given by a unique (up to sign) irreducible squarefree polynomial in Z[k, c]. We use the 
symbol A; nt also for that polynomial. To be precise, A jnt is the discriminant of a number 
held L with K D L D Q, namely L is the held of definition of the hnite ft-scheme Q. 

Next, for any i G {1,2,... ,19} consider the intersection of the steady state variety S with 
the hyperplane {ay = 0}. The Zariski closure of the image of S D {ay = 0} under the map 
7Tk, c is a hypersurface in C^ 9 x C 91 , dehned over Q, and we write A Xi=0 for the unique (up 
to sign) irreducible polynomial in Z[k, c] that vanishes on that hypersurface. We now define 

A . A; n t • 1cm (A a , 1= o , A.j; 2= o , ... , A 3 , ig= o ). 


This product with a least common multiple (1cm) is the discriminant for our problem. 

Example 4.3. The degree of A int as a polynomial only in c = ( 01 , 02 , 03 , 04 , 05 ) equals 34. 
To illustrate this, we set c = (5,16 + C, | — C, | + ( 7 , 3 — ( 7 ) where C is a parameter, and 


'9 922 4 22 44131 77 

k = ( —, —, 3, , 3,1,1,100, -, 80,100,1,3, —, —, 38, —, —,4, —, —, 1, —, 19, ——, 1,1,5,1 

5533 5 33 55852 44 


Under this specialization, the polynomial Ai nt becomes an irreducible polynomial of degree 
34 in the parameter C. Its coefficients are enormously large integers. It has 14 real roots. 

For the other factors A Xi= o of the discriminant, we find the following specializations: 


x\ — > 0, X 2 —> 0, X 3 —> 0, xy — > (C+16)(5(7—8), X 5 —> (7+16, x% -+ ((7+16)(5(7+6), 
xj —> (7+16, X 8 —> 5(7 — 8 , xg —> 5(7 + 6 , x\q — > a quartic q(C), x\\ —> 0 , X 12 —> (7—3, 
X13 -+ (7-3, X14 -> (<7+16)(5(7—8 ), x 15 -+ (<7+16)(5C+6), x 16 -+ (<7+16)(5(7—8 ) , 
X 17 -+ (<7+16)(5(7+6), xi 8 -+ ((7+16)(5(7—8 )q(C), xig -+ ((7+16)(5(7+6). 


These polynomials have 8 distinct real roots in total, so the total number of real roots of the 
discriminant is 14 + 8 = 22. These are the break points where real root behavior changes: 


(9,0) 

-77.2388 

(9,0) 

-16.0000 

(9,0) 

-5.28669 

(7,0) 

-1.57472 

(9,0) 

-1.46506 

(9,0) 

-1.34899 

(7,0) 

-1.29581 

(9,0) 

-1.20000 

(9.1) 

-1.19215 

(9,1) 

-1.18389 

(7,1) 

-0.584325 

(9,3) 

-0.361808 

(7.3) 

0.191039 

(5,1) 

1.30812 

(7,1) 

1.33197 

(5,1) 

1.60000 

(5,0) 

1.60161 

(3,0) 

3.0000 

(3,0) 

4.26306 

(5,0) 

11.1174 

(7.0) 

21.4165 

(9,0) 

310.141 

(9,0) 





In this table, we list all 22 roots of the specialized discriminant A(C). The eight boldface 
values of C are the roots of Q: here one of the coordinates of x becomes zero. At the other 
14 values of ( 7 , the number of real roots changes. Between any two roots we list the pair 
(r, p ), where r is the number of real roots and p is the number of positive real roots. For 
instance, for —0.361808 < C < 0.191039, there are 7 real roots of which 3 are positive. 
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5. Algebraic Matroids and Parametrizations 


Question 5 asks: Suppose we can measure only a subset of the species concentrations. Which 
subsets can lead to model rejection? This issue is important for the Wnt shuttle model 
because, in the laboratory, only some of the species are measurable by existing techniques. 


We shall address Question 5 using algebraic matroids. Matroid theory allows us to analyze 
the structure of relationships among the 19 species in Table [I] This first appeared in j IT . 
We here present an in-depth study of the matroids that govern the Wnt shuttle model. 

An introduction to (algebraic) matroids can be found in 21 ; they have been applied in 


14,15 to problems involving the completion of partial information. General algorithms for 
computing algebraic matroids are derived in 1241. We briefly review basic notions. 


Definition 5.1. A matroid is an ordered pair (A", X), where X is a finite set, here regarded as 
unknowns, and X is a subset of the power set of X. These satisfy certain independence axioms. 
For an algebraic matroid, we are given a prime ideal P in the polynomial ring K[X] generated 
by X, and X consists of subsets of X whose images in K[X]/P are algebraically independent 
over K. Thus, the collection of independent sets is X = \Y Cl : P n K\Y] = {0}}. 

1. Bases are maximal independent sets, i.e. subsets in X that have maximal cardinality. 


2. Rank is a function p from the power set of A" to the natural numbers, which takes as input 
a set Y C A" and returns the cardinality of the largest subset of Y in X. 

3. Closure is a function from the power set 2 X to itself. The input is a set Y and the output 
is the largest set containing Y with the same rank. 

4. Flats are the elements in 2 A that lie in the image of the closure map. 

5. Circuits are the sets of minimal cardinality not contained in X. 


We are here interested in the matroid that is defined by the prime ideal P = I m in Q(k) [x]. Its 
ground set X is the set of species concentrations {aq,..., £ 19 }. Since V(I m ) is 5-dimensional, 
each basis consists of five elements in X. In our application, bases are the maximal subsets of 
X that can be specified independently at steady state; they are also the minimal-cardinality 
sets that can be measured to learn all species concentrations. The rank of a set Y indicates 
the number of measurements required to learn the concentrations for every element of Y. 
Flats are the full subsets that are specified by any given collection of measurements. 

Circuits furnish our answer to Question 5: they are minimal sets of species that can be used 
to test compatibility of the data with the model. For each circuit Y there is a unique-up- 
to-scalars relation in I m fl Q(k)[F], called the circuit polynomial of Y. If the measurements 
indicate that this relation is not satisfied, then the model and data are not compatible. 

Proposition 5.2. The algebraic matroid of I m has rank 5. It has 951 circuits, summarized 
in Table [5| Of the 11628 subsets of X of size 5, precisely 2389 are bases. The 2092 bases 
summarized in Table\6\have base degree 1, while the remaining 297 have base degree 2. 


The computation of this matroid was carried out using the methods described in |24]. It was 
first reported in 17 , along with the matroids of alternative models for the Wnt pathway. 


The idea there was to find subsets of variables that were dependent for different models. 
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Our matroid analysis here goes beyond 17 in several ways: 


1. We keep track of the parameters k. We take our circuit polynomials to have (relatively 
prime) coefficients in Z[k]. This gives us a new tool for model rejection, e.g. in situations 
where only one data point is known but some parameter values are available. 


2. We show how circuits can be used in parameter estimation; this will be done in Section [8j 

3. We use the degree-1 bases to derive rational parametrizations of the variety V(I m ). 

We now explain Table [5] A circuit polynomial has type (■ i,j ) if it contains i species concen¬ 
trations (x-variables) and j rate parameters (k-variables). The entry in row i and column j 
in Table [5] is the number of circuits of type (i,j). Zero values are omitted for clarity. 



2 

3 

4 

5 

6 


2 3 4 

5 

6 


2 3 4 

5 

6 

2 

5 

1 




12 

13 

10 


22 

8 

58 


3 


6 




13 

13 

15 

2 

23 

4 

56 

5 

4 

1 

5 




14 

19 

16 

1 

24 


54 

14 

5 


6 

1 



15 

17 

21 

4 

25 


53 

15 

6 


7 

5 



16 

15 

11 

2 

26 


8 

16 

7 


5 

3 



17 

16 

32 

9 

27 

12 

56 

16 

8 


1 

11 

1 


18 

4 

6 

2 

28 

2 


2 

9 


6 

12 

3 


19 

26 

36 

11 

29 


29 

14 

10 



11 

1 


20 

44 

1 

1 

30 




11 


4 

7 

11 

1 

21 

26 

27 

9 

31 
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TABLE 5. The 951 circuit polynomials, by numbers of unknowns ay and kj. 


Example 5.3. There are five circuits of type (2, 2). One of them is ay = — kiXi + k 2 X 2 - Most 
of the 951 circuit polynomials in I m are more complicated. In particular, they are non-linear 
in both x and y. For instance, the unique circuit polynomial of type (6,11) equals 

( — ^ 15 ^ 17 ^ 19 ^ 20^25 — ^ 16 ^ 17 ^ 19 ^ 20 ^ 25 )^ 7 ^ 9^13 
+ (^14^16^18^21^24 + ^14^16^19^21^24 + ^14^16^18^22^24 + &14&16&19^22&24)x3Xl2ay9. 

In Section [7J we will consider the role of these nonlinear functions in parameter estimation. 

Given a basis Y of an algebraic matroid, its base degree is the length of the generic fiber of 
the projection of V(P) onto the E-coordinates (cf. 124]). Bases with degree 1 are desirable: 

Proposition 5.4. Let P C K[X] be a prime ideal, Y a basis of its algebraic matroid, 

|AG| = n, and |Y| = r. If Y has base degree 1 then V(P) is a rational variety, and the basic 
circuits ofY specify a birational map tpy '■ K r — ■> K n whose image is Zariski dense in V(P) 

Proof. For each coordinate ay in X\Y there exists a circuit containing Y U {ay}; this is the 
basic circuit of (Y, ay). Since Y has base degree 1, the generic fiber of the map V(P) —> K r 
consists of a unique point. Therefore the circuit polynomial is linear in ay. It has the form 

Pi(Y) ■ ay + qi(Y), where p h G K[Y]. 

The ^-coordinate of the rational map (fy equals ay if ay G Y and — qi(Y)/pi(Y ) if ay ^ Y. □ 
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From Propositions 


5.2 


and 


5.4 

kF- 


we obtain 2092 rational parametrizations of the variety V(I m ). 
K 19 , where Y runs over all bases of base degree 1. Using 


These are the maps p y 
these <py, we obtain 2092 representations of the steady state variety Q as a subset of K~\ 
where now K = Q(k, c). Namely, we consider the preimages of the five hyperplanes defined 
by <§■ These are hypersurfaces in K 5 whose intersection represents the nine points in Q. 
We performed the following computation for all 2092 bases Y = {y 4 ,..., y 5 } of base degree 1: 

1. Substitute x = ipyiv i, •••, 1/5) into the hve linear equations (|2]) . 


2. Clear the denominators d 4 , ..., d$ in each equation to get polynomials hi,... ,h 5 in Y. 

3. The saturation ideal Jy = {hi, ..., /15) : {did 2 • • ■ d§)°° represents the preimage of Q. 


Given such a wealth of parametrizations, we seek one where Jy has desirable properties. 
We use the following criterion: consider subsets of hve of the generators of Jy, compute the 
mixed volume of their Newton polytopes, and fix a subset minimizing that mixed volume. In 
the census of 2092 bases in Table [6j that minimum is referred to as the mixed volume of Y. 


Mixed Volume 

5 

9 

10 

11 

12 

13 

14 

15 

16 

20 

23 

24 

25 

30 

35 

42 

45 

Frequency 

2 

416 

6 

73 

50 

167 

563 

751 

10 

12 

6 

1 

11 

12 

4 

4 

4 


Table 6. Reducing the steady state equations to the 2092 bases of base degree 1 


By Bernstein’s Theorem, the mixed volume is the number of solutions to a generic system 
with the hve given Newton polytopes. We seek bases Y where this matches the number nine 
from Theorem 11.11 We see that the mixed volume is nine for 416 of the bases in Table [HI 

Example 5.5. The basis Y = {xi,x±,xq,x&,xiz} has base degree 1 and mixed volume 9. 
The remaining variables can be expressed in terms of Y as follows. For brevity, we set 

r{x 4 ,x 6 ) = k 9 kuk 20 k 22 x 4 x 6 + k 9 k u {k 21 T k 22 )(k 23 + k 31 )x 4 

+ ^20^22(^10 + ku)(kl 3 + k3o)x6 + (fcio + kii){k-2l + ^22X^13^23 + ^23^30 + & 13 & 3 l)- 


X 3 

x 5 

x 7 

Xg 

Xn 


kik 26 

r“7— 

k 2^27 

k 4 k 3 k 5 (k 7 + k 8 ) x 4 x 4 
k 2 k 6 k 8 (k 4 + k 5 ) x 8 
kik 3 k 5 k 2S {k 7 + fc 8 ) 3:1X4 
k 2 kgk 8 k 2 g {k 4 + k 5 ) x 8 

fc 6 fc 8 fci 4 fci 6 fc26fc29(fc4+A-5)(fci8+fcl9) 
kzk^k^k\’jk\Qk2-rk2s{.^+kp) (fci.5+^16) 

fcl2(fcl0+fcll)(fc2Qfc22^6 + (fc21+fc22K^23+fc3l)) 

r(x 4,Xe) 


Xi2 

Xu 

X 15 

Xl 6 

X U 

Xl 8 


^12^30(^10 + kn)(k 2 i + k 22 ) 
r(x 4 ,x 6 ) 


X 19 


_ r(x4,xe) _ fc 25 

fcl2fe3o(felO+fcll)(fc21+fe22) 13 

kik 3 

I 71 A 1 \ XlXi 

k 2 [k 4 + k 5 ) 
kikuk 2 6 
k- 2 k 27 {ki 5 + kio) 
kik 3 k 5 

k 2 k s {k 4 + k 5 ) XlX4 
k\k\ 4 k\^k 2 Q 
k 2 kigk 27 {ki 5 + kio) 

kgki2(k2ok22X(i + (k2l+k22)(k2S+k3l)) 

r(x4,x 6 ) 4 


^12^20^30(^10 + ku) 
r(x 4 ,x 6 ) 


Xg 
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This map ipy is substituted into (j2]), and then we saturate. The resulting ideal Jy equals 
{ayXQXg + a 2X4 + CX3XQ, CX4X1XQ + CX5X1 + ol gXg + ay, 

(23X1X4 + a gXs + quo , (211X4X6X13 + 01.12X4X13 + (213X^X13 + (214X13 + (215, 

(2iqX4Xq + (2ijXq + 0(18^4^6+ a 19 x q + 020^8 A a 2lXl + 0^22^4 + 0^23^6 + ®24^8 + Ol2§) i 

where the a 1 ,... ,025 are certain explicit rational functions in the k-parameters. 


6 . Polyhedral Geometry 


Dynamics of the system while not at steady state cannot typically be studied with algebraic 
methods. One exception is the set of all possible states accessible from a given set of initial 
values via the chemical reactions in the model. This set is called a stoichiometric compatibility 
class in the biochemistry literature. Mathematically, these classes are convex polyhedra. We 
determine them all for the Wnt shuttle model. This resolves Problem 6 from the Introduction. 

The conservation relations ([ 2 J) define a linear map x from the orthant of concentrations M> 9 0 
to the orthant of conserved quantities M> 0 - We express this projection as a 5 x 19-matrix: 


( 7 ) 



1 1 
1 1 


. . . A 

1111 
1 . . . 

. 1 . . 

■ ■ ■ 7 


f xi \ 

X2 

X3 

XlS 

\X19J 


Let P c denote the fiber of the map x f° r c *= ^>o- This is known in the biochemical 
literature as the invariant polyhedron or the stoichiometric compatibility class of the given x; 
see e.g. 25, (3)]. The fiber over the origin is Po = BGojeio, e n}, the two-dimensional 
orthant formed by all positive linear combinations of eio and en. If c G M> 0 is an interior 
point, then P c is a 14-dimensional convex polyhedron of the form Po x P c where P c is 
a 12-dimensional (compact) polytope. Two vectors c and c' are considered equivalent if 
their invariant polyhedra P c and P c / have the same normal fan. This property is much 
stronger than being combinatorially isomorphic. The equivalence classes are relatively open 
polyhedral cones, and they define a partition of M>o- This partition is the chamber complex 
of the matrix (J?J). For a low-dimensional illustration, see 125, Figure 1], Informally speaking, 
the chamber complex classifies the possible boundary behaviors of our dynamical system. 


Proposition 6.1. The chamber complex of our 5x19 -matrix divides M> 0 into 19 maximal 
cones. It is the product of a ray, M>o, and the cone over a subdivision of the tetrahedron. 
That subdivision consists of 18 smaller tetrahedra and 1 bipyramid, described in detail below. 


Proof. The product structure arises because the matrix has two blocks after permuting 
columns, an upper left 4x17 block and a lower right 1x2 block (1 1). Our task is to compute 
the chamber decomposition of R> 0 defined by the 4 x 17-block. After deleting zero columns 
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and multiple columns, we are left with a 4 x 7-matrix, given by the seven left columns in 


M = 


a b c d e f g h i j k l 

/0 1 0 0 1 0 0 0 1 1 1 l\ 

1111000 11112 
0010010 10111 
yo ooiooi i i o i \) 


The correspondence between the seven left columns of M and the columns of (Jt]) is as follows: 

a = {X4, £5, Xq, X 7 , £18, £19}, 6 = {^14; ^15}) c = {x 16}, 

d = {x 17 }, e = {xi, x 2 , 2:3}, / = K}, 9 = {^ 9 }- 

The remaining columns of M are additional vertices in the subdivision. 

The following table lists the 19 maximal chambers. For each chamber we list the extreme rays 
and the facet-defining inequalities. For instance, the chamber in M > 0 denoted by efjk is the 
orthant spanned by the columns e, /, j and k of the matrix M times the ray (0, 0,0, 0, 1) T . It 
is defined by C 5 > 0 together with the four listed inequalities: C 4 > 0, min(ci, C 3 ) > C 2 > C 4 . 


abed 

bcdl 

efgk 

bcjl 

bdil 

beij 

cdhl 

cfhj 

dghi 

egik 

fghk 

efjk 

bijl 

chjl 

dhil 

ghik 

eijk 

fhjk 

hijkl 


{c 2 


{C4, C3, Ci, C 2 C4 C3 Ci} 

{c 2 — C3 — Cl, C 2 — C4 — Cl, C 2 — C4 — C3, —C 2 + C4 + C3 + Cl} 

(c 2 , —C2 + C4, —C2 + C3, —C2 + Cl} 

(c4, —C2 + Cl + C3, C2 — C3 — C4, C2 — Cl — C4} 

{ c 3; ~ c 2 + C4 + Ci, C2 — C4 — C3, C2 — C3 — Ci} 

{C3, C4, Ci — C 2 , C2 — C3 — C4} 

{ci, —C2 + C3 + C4, C2 — C4 — C3, C2 — C4 — Ci} 

{C4, Ci, —C 2 + C3, C2 — C4 — Ci} 

{Cl, C3, —C 2 + C4, C2 — Ci — C3} 

(C3, —C2 + C4, C 2 — C3, —C2 + Ci} 

{Ci, —Ci + C2, —C2 + C3, —C2 + C4} 

{C4, Ci — C2, C 2 — C4, —C2 + C3} 

{c 2 — Ci, —C2 + Ci + C3, —C2 + C4 + Ci, C2 — €3 — C4} 

{c 2 — C3, —C2 + C4 + C3, —C2 + C3 + Ci, C2 — C4 — Ci} 

{c 2 — C4, —C2 + C4 + Cl, —C2 + C3 + C4, C2 — Cl — C3} 

(C4 — C2, C2 — C3, C2 — Ci, —C2 + Ci + C3} 

(c 2 — C4, C2 — C3, Ci — C 2 , —C2 + C3 + C4} 

(c 2 — C4, C2 — Cl, —C2 + C3, —C2 + C4 + Cl} 

- C4, C2 — C3, C2 — Cl, —C2 + C4 + C3, —C2 + C4 + Cl, —C2 + C3 + Cl} 


Interpreting the columns of M as homogeneous coordinates, the table describes a subdivision 
of the standard tetrahedron into 18 tetrahedra and one bipyramid hijkl. These cells use the 
12 vertices a,b,... ,1. The reader is invited to check that this subdivision has precisely 39 
edges and 47 triangles, so the Euler characteristic is correct: 12 — 39 + 47 — 19 = 1. □ 


We shall prove the following result about the Wnt shuttle model. 


Proposition 6.2. Suppose that the rate constants ki and the conserved quantities Cj are all 
strictly positive. Then no steady states exist on the boundary of the invariant polyhedron P c . 
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Proof. Consider the two components I m and I e of the steady state ideal I given in Lemma 3.1 


We intersect each of the two varieties with the affine-linear space defined by the conservation 
relations (J2|) for some c G R> 0 - We claim that all solutions x satisfy x* ^ 0 for i — 1, 2,..., 19. 

For the main component V(I m ), we prove this assertion with the help of the parametrization 
ip Y from Example 5.5 If the values of x 4 , x 4 , x 6 , x 8 , x 13 and of the expression r(x 4 ,x 6 ) are 


nonzero, then each coordinate of ip Y is nonzero. We next observe that r(x 4 ,x 6 ) > 0 for any 
k > 0 and x > 0 . A case analysis, using binomial relations in the ideal J m , reveals that if 
any of x 4 , x 4 , Xq, x 8 , X13 are zero, some coordinate of c is forced to zero as well: 


Xl 

= 0 


X 2 , X 3 , X 44 , X 15 

= 0 


Cl 

= 0 

Xl3 

= 0 


X 12 

= 0 


c 5 

= 0 

x 4 

= 0 

=> 

X 5 , X 6 , X 7 , Xi 4 ,Xi5,Xi6,Xi7,Xi 8 ,Xi9 

= 0 


c 2 

= 0 



or 

XS: X\q 

= 0 


c 3 

= 0 

x 6 

= 0 


X 9 , X 17 

= 0 


c 4 

= 0 



or 

X 4 

= 0 


c 2 or c 3 

= 0 

X8 

= 0 

=> 

Xl6 

= 0 


c 3 

= 0 


It remains to consider the extinction component. Its ideal I e contains the set b U l — 
{x 4 , X 2 , X 3 , x 44 , X 15 }. The corresponding columns of the matrix in ([7]) are the only columns 
with a nonzero entry in the fourth row. This implies that c 4 = 0 holds for every steady state in 
V (/ e ). We conclude that there are no steady states on the boundary of the polyhedron P c . □ 

Remark 6.3. In this proof we did not need the detailed description of the chamber complex, 
because of the special combinatorial structure in the Wnt shuttle model. In general, when 
studying chemical reaction networks that arise in systems biology, an analysis like Proposition 
| 6 . 1 | is requisite for gaining information about possible zero coordinates in the steady states. 


7. Parameter Estimation 


Question 7 asks: What information does species concentration data give us for parameter 
estimation? This question is of particular importance to experimentalists, as species concen¬ 
trations depend on initial conditions, whereas parameter values are intrinsic to the biological 
process being modeled. Identihability of parameters has been studied in many contexts, no¬ 
tably in statistics ( 8 ) and in biological modeling 19 . Sometimes, as in fl9), parameters are 


determined from complete time-course data of the dynamical system, making a differential 
algebra approach desirable. In the present paper we focus on the steady state variety, so we 
consider data collection only at steady state. We assume that there is a true but unknown 
parameter vector k* G R 31 of rate constants, and our data are sampled from the positive 
real points x on the variety in R 19 that is defined by the 19 polynomials in ([!]). 


7.1. Complete Species Information. The first algebraic question we answer: To what 
extent is the true parameter vector k* determined by points on its steady state variety? 

To address this question, we form the polynomial matrix F(x) of format 19 x 31 whose 
entries are the coefficients of the right-hand sides of ([!]), regarded as linear forms in k. With 
this notation, our dynamical system ([TJ) can be written in matrix-vector product form as 

= F(x) • k. 
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Our data points are sampled from 

(8) { x G !> 0 : F(x) • k* = 0 }. 

Let xi,X 2 ,X 3 ,... denote generic data points in (J8]) . The set of all parameter vectors k that 
are compatible with these data is a linear subspace of M 31 , namely it is the intersection 

(9) kernel(F(x 1 )) D kernel(F(x 2 )) fl kernel (F(x 3 )) D • • • 


The best we can hope to recover from sampling data is the following subspace containing k*: 


( 10 ) 


P'j kernel(F(x)) C M 31 . 

x in |8| 


We refer to (10) as the space of parameters compatible with k*. A direct computation reveals: 


Proposition 7.1. The space of all parameters compatible with k* is a 14 -dimensional sub¬ 
space o/M 31 . //x is generic then the kernel of Ffx) is a 17 -dimensional subspace ofM. 31 . 


This has the following noteworthy consequence for our biological application: 

Corollary 7.2. The parameters of the Wnt shuttle model are not identifiable from steady 
state data, but there are 14 degrees of freedom in recovering the true parameter vector k*. 


Our next step is to gain a more precise understanding of the subspaces in Proposition 7.1 


To do this, we shall return to the combinatorial setting of matroid theory. We introduce two 
matroids on the 31 reactions in Table [2] The common ground set is Ii = {k\, k 2 , • • •, fc 3 i}. 
The one-point matroid A4 one is the rank 17 matroid on K defined by the linear subspace 
kernel (F(x)) of M 31 where x e M 19 is generic. The parameter matroid Al par is the rank 14 
matroid on K defined by the space (10) of all parameters compatible with a generic k*. The 
following result, obtained by calculations, reflects the block structure of the matrix F(x). 


a) Single Data Point b) Multiple Data Points 



FIGURE 1. Graphic representation of the one-point matroid M. one of rank 17. 

The rank 4 component of the rank 14 parameter matroid A4 par is not graphic. 

Proposition 7.3. The one-point matroid A4 one is the graphic matroid of the graph shown in 
Figure^ 7] a). Its seven connected components are matroids of ranks 3, 3, 7,1,1,1,1. The rank 
14 parameter matroid A4 par is obtained from A4 one by specializing the rank 7 component to 
the rank 4 matroid on 11 elements whose affine representation is shown in Figure^ 7] b). 
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This characterizes the combinatorial constraints imposed on the parameters k by measuring 
the species concentrations at steady state. For a single measurement x, the result on A4 one 
tells us that the 19 x 31-matrix F{x) has rank 14 = 31 — rank(A4 one ). After row operations, 
it block-decomposes into two matrices of format 3x6, one matrix of format 4 x 11, and four 
matrices of format 1x2. Each of these seven matrices is row-equivalent to the node-edge 
cycle matrix of a directed graph, with underlying undirected graph as in Figure [I] (a). 

Consider the graph with edges 9,10,11,12,13, 20, 21, 22,23, 30, 31. The cycle {22, 23, 30, 31} 
reveals that our measurement x imposes one linear constraint on k 2 2 , ^ 23 , & 30 , ^31 • If we take 
further measurements, as in (J9]) , then six of the seven blocks of F(x) remain unchanged. Only 
the 4 x 11-block of F(x) must be enlarged, to a 7 x 11-matrix. The rows of that new matrix 
specify the affine-linear dependencies among 11 points in M 3 . That point configuration 
is depicted in Figure [I] (b). For instance, the points {9,10,11} are collinear, the points 
{20, 21, 22} are collinear, but these two lines are skew in M 3 . From the other line we see that 
that repeated measurements at steady state impose two linear constraints on k 2 2 , k 2 3 , k 30 , Am- 


7.2. Circuit Data. The second question we address in this section: Given partial species 
concentration data, is any information about parameters available? In Section 7.1, all 19 
concentrations Xi were available for a steady state. In what follows, we suppose that Xi can 
only be measured for indices i in a subset of the species, say C C {1,..., 19}. In our analysis, 
it will be useful to take advantage of the rank 5 algebraic matroid in Proposition |5.2[ since 


that matroid governs dependencies among the coordinates aq,..., aqg at steady states. 


We here focus on the special case when C is one of the 951 circuits of the algebraic matroid 
of I m . Let fc be the corresponding circuit polynomial, as in Table [ 5 } We regard fc as a 
polynomial in x whose coefficients are polynomials in Q[k]. Suppose that fc has r monomials 
x ai ,..., x“ r . We write F c e Q[k] r for the vector of coefficients, so our circuit polynomial is 
the dot product fc( k, x) = Fc-(k) • (x Ql ,..., x ar ). We write Vc C M r for the algebraic variety 
parametrized by F c (k). Thus Vc is the Zariski closure in M r of the set {F c (k') : k' G M 31 }. 

Our idea for parameter recovery is this: rather than looking for k compatible with the true 
parameter k*, we seek a point y = Fc(k) in Vc that is compatible with F c (k*). And, only 
later do we compute a preimage of y under the map M 31 —> M r given by Fc- Most interesting 
is the case when Vc is a proper subvariety of W. Direct computations yield the following: 


Proposition 7.4. For precisely 288 of the 951 circuits C of the algebraic matroid of the 
steady state ideal I m , the coefficient variety Vc is a proper subvariety in its ambient space 
R r . In each of these cases, the defining ideal ofVc is of one of the following four types: 


(11) (2/22/6 - 2/32/s) 

(12) (2/52/6 - 22/32/7,2/1 - 42/22/7,2/32/5 - 2y 2 y G , 2/22/1 - 2/32/7) 

(13) (2/32/5 - 2/22/52/6 + 2/12/1) 

(14) ( 22 / 3 2/4 - 2 / 22 / 5 , 2 / 22/3 - 2//i2/5, y\ ~ tyilu) 
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Example 7.5. Consider the circuit C = {6,10,18}. The circuit polynomial fc equals 

(^13^20^22 + ^20^22^30) ' ^6^10 + ^11^20^22 ' ^6^18 ~~ ^12^20^22 ' ^6 
+ (^13^21^23 + ^13^22^23 + ^21^23^30 + ^22^23^30 + ^13^21^31 + ^13^22^3l) ' ^10 
+ (^11^21^23 + ^11^22^23 + &11&21&31 + ^11^22^3l) ' ^18 
— (^12^21 ^23 + ^12^22^23 + ^12^21^31 + ^12^22^3l)- 

Here r = 6 and we write -Fc(k) = ( 2 / 1 , 2 / 2 , 2 / 3 , 2 / 4 , 2 / 5 ^ 2/6) for the vector of coefficient polynomi¬ 
als. The variety Vc is the hypersurface in M 5 defined by the equation 2/22/6 = 2 / 32 / 5 . 

We now sample data points x,; from the model with the true (but unknown) parameter vector 
k*. Each such point defines a hyperplane (y 6f r : y • (x® 1 ,... , x® r ) = 0}. The parameter 
estimation problem is to find the intersection of these data hyperplanes with the variety Vc- 
That intersection contains the point y* = Ec(k*), which is what we now aim to recover. 


7.3. Noisy Circuit Data. The final question we consider in this section is: Given partial 
species concentration data with noise, is any information about parameters available? 

As in Section 7.2, we fix a circuit C of the algebraic nratroid in Section [5j and we assume 
that we can only measure the concentrations x 3 where j G C. Each measurement x* G M c 
still defines a hyperplane y ■ (x® 1 ,... , x“ r ) = 0 in the space M r . But now the true vector 
y* = Ec(k*) is not exactly on that hyperplane, but only close to it. Hence, if we take s 
repeated measurements, with s > r, the intersection of these hyperplanes should be empty. 

We propose to find the best fit by solving the following least squares optimization problem: 


(15) Minimize ^(y ■ (x® 1 ,... , x“ r ) ) 2 subject to y G Vc fl § r \ 

2=1 


where § r_1 = (y G M r : y\ + y\ + • • • + y 2 r = 1} denotes the unit sphere. When the variety 
Vc is the full ambient space M r , this is a familiar regression problem, namely, to find the 
hyperplane through the origin that best approximates s given points in M r . Here “best” 
means that the sum of the squared distances of the s points to the hyperplane is minimized. 
This happens for 663 of the 951 circuits C, and in that case we can apply standard techniques. 


However, for the 288 circuits C identified in Proposition |7.4[ the problem is more interesting. 
Here the hyperplanes under consideration are constrained to live in a proper subvariety. I 11 
that case we need some algebraic geometry to reliably find the global optimum in (15). 

Our problem is to minimize a quadratic function over the real affine variety Vc 0 S r_1 . The 
quadratic objective function is generic because the x* are sampled with noise. The intrinsic 
algebraic complexity of our optimization problem was studied by Draisma et al. in (5]. That 
complexity measure is the ED degree of Vc 0 S”" 1 , which is the number of solutions in C r 
to the critical equations of (15). Here, by ED degree we mean the ED degree of Vc n§ r_1 , 
when considered in generic coordinates. This was called the generic ED degree in 1201. 


We illustrate our algebraic approach by working out the first instance 0) in Proposition |7.4| 


Example 7.6. Suppose we are given s noisy measurements of the concentrations xe, x\o, xis- 
In order to find the best fit for the parameters k, we employ the circuit polynomial fc in 
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Example 7.5 We compute y G M b by solving the corresponding optimization problem ( 16 ). 


This problem is to minimize a random quadratic form subject to two quadratic constraints 
( 16 ) 


« - </ 3 </S 


2 , 2 , 2 , 2 , 2 , 2 
— Di + I/2 + 2/3 + 2/4 + Vs + % 


1 = 0 . 


We solve this problem using the method of Lagrange multipliers. This leads to a system of 


polynomial equations in y. Using saturation, we remove the singular locus of (16), which is 
the circle {y G M 6 : y\ + y\ — 1 = 2/2 = 2/3 = 2/5 = 2/6 = 0}. The resulting ideal has precisely 
40 zeros in C 6 . In the language of (5,20 , the generic ED degree of the variety (16) equals 40. 


8 . From Algebra to Biology 

The aims of this paper are: (1) to demonstrate how biology can lead to interesting ques¬ 
tions in algebraic geometry, and (2) to apply new techniques from computational algebra 
in biology. So far, our tour through (numerical) algebraic geometry, polyhedral geometry 
and combinatorics has demonstrated the range of mathematical questions to explore. In 
this section, we will focus on translating our analysis into applicable considerations for the 
research cycle in systems biology, which is illustrated in Figure [2j In what follows we discuss 
some concrete applications and results pertaining to the steps (a), (b) and (c) in Figure [ 2 } 


Biological 
hypothesis 
Wnt Signal 



Response? 


Model 

construction 

x i + x j x i + x i 

Model f(x,k) 
Steady states f-0 


Biological 

insight 


(a) Model analysis 


Theorem I 


im 


■Re 


Discriminant, Ex. 4.3 


Y (b) Experimental design 



Prop. 6.2, Ex. 8.1 


Matroids 
Prop. 5.2, Ex. 8.2 


Parameter estimation, Cor. 7.2, Ex. 8.4 

Model rejection via Circuits Model variety 
Ex 5.3, Ex. 8.3 ^--- '• 


*) 


Data point 


(c) Model and data compatibility 


Experiments 


Figure 2 . Systems biology cycle informed by algebraic geometry and com¬ 
binatorics. (a) Model analysis. See Sections 1, 3, 4. (b) Experimental design. 
See Sections 5 and 6. (c) Model and data compatibility. See Sections 5 and 7. 


Analysis of the Model: Before any experiments are performed, our techniques inform the 
modeler of the global steady-state properties of the model. The number of real solutions to 
system 0-0 , stated in Theorem 0 governs the number of observable steady states. Var¬ 
ious sampling schemes demonstrated that most parameter values lead to only one observable 
steady state. We produced a set of parameter values and conserved quantities with three 
real solutions, and two solutions are also attainable. If the “true” parameters k* and c* 
admit multiple real solutions, then multistationarity of the system is theoretically possible. 

If multiple states are observed experimentally, then the model must be capable of multista¬ 
tionarity. In the Wnt shuttle model, the system is capable of multiple steady states; however, 
based on parameter sampling, the frequency of this occurrence is low, and parameters in this 
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regime are somewhat stable under perturbation. The discriminant of the system is a poly¬ 
nomial of degree 34 in c, and our analysis along a single line in c-space illustrates the high 
degree of complexity inherent in the full stratification of the 36-dimensional parameter space. 

Experimental Design: In Section[ 6 j the combinatorial structure of the various stoichiomet¬ 
ric compatibility classes was fully characterized. As the conserved quantities c = (ci,..., c 5 ) 
range over all positive real values, the set of all compatible species-concentration vectors x 
will take one of 19 polyhedral shapes P c . This may find application in identifying multiple 
steady state solutions for specific rate constants k. A natural choice for initial conditions 
when performing experiments is on or near the vertices of the 14-dimensional polyhedron P c . 

Example 8.1. Suppose the conserved quantities vector lies in the bipyramid, e.g. c = 
(1,2, 2, 2,3). The preimage of c in x-space is a product of the orthant M>o{ei 0 ,en} and 
a 12-dimensional polytope with 400 vertices: (1, 0, 0, 2, 0, 0, 0, 2, 2, 0,0, 3, 0, 0,0, 0, 0, 0,0, 0), 
and 399 of its permutations. This product is the polyhedron P c . If we have control over initial 
conditions, beginning near the vertices positions us to find interesting systems behavior. 


In the laboratory, the experimentalist makes choices of what to measure and what not to 
measure. For instance, measuring a particular Xj may be infeasible, or there may be a 
situation in which measuring concentration Xj can preclude measuring concentration Xj. 

For every strategy, we fix a cost vector , listing the costs of making each measurement. We use 
the symbol N to indicate infeasible measurements. Suppose there are two different ways to 
run the experiment; then we have a 2 x 19 cost matrix P , whose rows are cost vectors for each 
experiment. We multiply P by the 0-1-incidence matrix for the 951 circuits of Proposition 
5.2| That matrix has a 1 in row i and column j if circuit j contains species i, and 0 otherwise. 
The product is a matrix of size 2 x 951. For N —> oo, the 2 x 951 matrix has a finite entry 
in position (i,j) precisely when the strategy i can measure the circuit j. Minimizing over 
those finite cost entries selects the most cost-effective experiment to measure a circuit. 


Example 8.2. Suppose that none of the intermediate complexes x\ 3 ,..., x 1 9 are measurable, 
and that we are able to measure only one Phosphatase concentration (x 4 or xs) in each 
experimental setup. A corresponding cost matrix might look like 


P 


1 1 1 N 1 1 1 1 1 1 1 1 iV iV iV iV iV iV iV 

111 1 llliVlllliViViViViViViV 


Multiplying by the circuit support matrix of size 19 x 951 reveals 82 feasible experiments: 
50 using the first row of P, and 32 using the second. With more refined cost assignment, 
this would decide not only feasibility but also optimal cost. In this way, the matroid allows 
us to choose cost-minimal experiments to obtain meaningful information for the model. 


Model and data compatibility: After an experiment is performed, the task of the modeler 
is to test the data with the model. One possible outcome is model rejection. If the data are 
compatible, then another outcome is parameter estimation. Both may provide insights for 
biology. The role of algebraic geometry is seen in (9,10 and shown in the next two examples. 


Example 8.3 (Model Rejection). Suppose that rate parameters ki are all known to be 1, 
and that we have collected data for variables x 4 , x 4 , X 14 . The circuit polynomial is kik 3 xix^ + 
(—— k 2 k^)xi 4 , which specializes to x 4 x 4 — 2x 44 . If the evaluation of the positive quantity 
|x 4 x 4 — 2 x 44 1 lies above a threshold e, then we can reject the model as not matching the data. 
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Every circuit polynomial of the matroid is a steady state invariant ; depending on which 
experiment was performed, the collection of measured variables must contain some circuit. 
Even if one can measure all 19 species at steady state, it is not possible to recover all 31 kinetic 
rate constants, but we do have relationships that must be satisfied among parameters [l 6 |. 

Example 8.4 (Parameter Estimation). Suppose that rate parameters are unknown, and that 
we have collected data for x 6 , Xi 0 , aq 8 . The corresponding circuit polynomial f c is shown in 
Example |7.5[ We know that the coefficients of fc satisfy the constraint y 2 y 6 = 2 / 32 / 5 . Suppose 
our experiments lead to the following ten measurements for the vector (xq, £ 10 , Xi 8 ): 


{(.715335, 4.06778,14.6806), (.390982,4.83152,6.08251), (.706539,4.98107,3.83617), 

(.14316,4.30851,12.5809), (.995583,4.01222,15), (.413817,4.08114,14.902), (.232206, 3.38274, 23.3162), 
(.219045, 5.06008,3.67175), (.704106, 3.52804, 21.1037), (.648732, 3.6505,19.7008)} 


The data lead us to the following function to optimize in (15): 


57.2345 y{ + 376.181 2/12/2 + 801.672^ - 27.5625y 1 y 3 - 96.U29y 2 y 3 
+3.36521 y\ + 179 A9ym + 564.0342/ 2 2 / 4 - A2.729y 3 y 4 + 178.839 y\ + 564.034y 1 ?/5 
+2424.317/22/5 - 144.77/32/5 + 1054.492/42/5 + 2263.22/| - 42.729 t/i2/6 
-144.72/22/6 + 10.33 92/32/6 - 83.80 7 22/ 4 2 / 6 - 269.74 92+2/6 + IO 2 t\ 


The global minimum of this quadratic form on the codimension 2 variety (16) has coordinates 

2/i = 0.183472, 2/2 = 0.152416, y 3 = 0.959232, 2/4 = 0.038042, 2/5 = 0.00335267, 2/6 = 0.211. 

Given these values, one now has three degrees of freedom in estimating the nine parameters 
ki that appear in the circuit polynomial fc- The other ten coordinates of k are unspecified. 
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